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I.  INTRODUCTION 

In  this  paper  we  show  how  molecular  dynamics  can  be  used  in  a  simple  manner  to  com¬ 
pute  the  contours  of  electronic  absorption  bands.  The  motivation  for  this  work  is  to  compute 
transient  electronic  spectra  for  many  atom  systems.1- 2  We  could  equally  use  Monte  Carlo  or 
explicit  integration  over  coordinates  to  compute  equilibrium  electronic  absorption  bands.  How¬ 
ever,  molecular  dynamics  provides  the  only  means  to  calculate  these  band  contours  in  a  non¬ 
equilibrium  system.  Thus,  by  demonstrating  that  molecular  dynamics  yields  the  correct  absorp¬ 
tion  band  contour  for  the  equilibrium  case,  we  can  justify  using  this  technique  in  the  nonequili¬ 
brium  case.  We  have  already  elsewhere1-2  applied  the  methods  derived  in  this  paper  to  the  cal¬ 
culation  of  the  transient  electronic  absorption  spectra  from  a  chemical  reaction  occurring  in 
solution.  In  another  paper3  the  electronic  absorption  spectrum  for  an  initially  thermal  equili¬ 
brium  system  is  computed  by  a  time  dependent  wave  packet  technique,  and  the  result  for  the 
same  I2  molecular  test  case  can  be  compared  to  the  more  classical  technique  developed  here. 

In  Section  II  we  discuss  the  theory  which  we  will  use,  as  well  as  a  harmonic  quantum 
correction  appropriate  for  thermal  equilibrium.  Section  III  contains  a  comparison  with  other 
semidassical  approaches  for  electronic  spectra.  The  potential  energy  and  electronic  transition 
dipole  versus  intemuclear  distance  functions  used  to  compute  the  sample  l2  visible  absorption 
band  contours  are  presented  in  Section  IV.  Section  V  discusses  the  computational  method  and 
Section  VI  compares  calculated  and  measured  equilibrium  I2  gas  phase  band  contours.  Finally, 
in  section  VII  we  discuss  the  results  and  their  significance. 

II.  THEORY 

A.  Absorption  cross  section 

The  absorption  cross  section  cr(w)  can  be  expressed  in  the  usual  first  order  time  depen¬ 
dent  perturbation  theory  approach  as4-6 

<t(oj)  —  r? ~  (2.1) 

7lcn  ,./ 

in  which  &>  is  the  angular  frequency.  D  is  Planck's  constant  divided  by  2ir.  c  is  the  speed  of 
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light,  n  is  the  index  of  refraction,  p,  and  pt  are  the  probabilities  in  the  ensemble  of  initial  and 
final  states,  (f  I  is  a  final  state,  |/)  an  initial  state,  «  is  a  unit  vector  in  the  direction  of  the  elec¬ 
tric  field  of  tne  probe  light,  p  is  the  electric  dipole  operator 

M-1 9,t,  (2-2) 

J 

in  which  qs  is  the  charge  and  r;  the  position  of  the  j  th  particle,  and  the  argument  to  the  Dirac 
8  function  is  in  terms  of  the  angular  frequency  difference 

«/,  -  ( £/-£ )/*  (2-3) 

between  the  energies  £/  and  £  of  the  final  and  initial  states. 

For  an  equilibrium  system 

Pi-p/  -  P/(l  -  e~s*')  (2.4) 

where  0  =  (kBT)~l  in  which  kB  is  Boltzmann’s  constant  and  fthe  temperature.  We  use  the 
usual  linear  response4-*  conversion  to  the  Heisenberg  picture  and  closure  to  write 

<r(w)  -  (2ir)-'X.“d/  e"« '%>,  (r|[a-»(0)l(a*»(r)]|i)  (2.5) 

in  which  m(0)  and  ji(t)  are  the  Heisenberg  time  dependent  dipole  moment  operators  evaluated 
at  times  0  and  i  If  we  adopt  a  classical  picture  for  the  dipole  moment  p(r),  we  may  rewrite6 
Eq.  (2.5)  in  terms  of  the  spectral  density  operator  D  as 

<rU)  - (  *>[i  *m)  )  (2.6) 

in  which  ( )  indicates  an  ensemble  average  over  initial  states  (a  thermal  average  for  the  equili¬ 
brium  case)  and  the  spectral  density  operator  per  unit  angular  frequency  is  defined  as 

D[f]  3  D[fJ]  -  (2ir)"'  lim  / Jt  *-•*/(/)  J*_V  e‘“'/V)  (2.7a) 

-  (2ir)'1  Urn |  fir*  e""'/(f)  | J-  (2.7b) 

For  a  nonequilibrium  system  there  is  no  a  priori  relationship  between  p,  and  pf.  We  will 
assume  here  for  non-equilibrium  systems  that  pf  =  0.  This  means,  for  example,  that  we  may 
start  with  a  small  density  of  final  states  and  illuminate  weakly  enough  so  that  the  final  states 
aren’t  pumped  appreciably.  The  appropriate  linear  response  relationship  is  then 

(D(5-m1)  (2.8) 

Hen  '  ' 

in  which  ' '  now  indicates  an  average  over  a  nonequilibrium  ensemble. 

B.  Semidassicai  time  dependent  dipole  moment 

Mukamel,  Stem,  and  Warshel7  ®  discuss  a  variety  of  semidassicai  approaches  to  electronic 
spectra,  and  what  we  illustrate  here  is  perhaps  the  simplest  and  most  classical  of  the  possible 
choices.  The  semidassicai  picture  for  the  dipole  moment  which  we  use  takes  the  nuclear 
motion  to  be  classical  and  the  electronic  transition  as  a  harmonic  osdllator  whose  angular  fre¬ 
quency  n /,  is  given  by  the  energy  difference  between  initial  and  final  electronic  state  potential 
energy  curves 

n/(-r'  lK,(rv)-  K((rv)l.  x  (2.9) 

in  which  Vf  (rv)  is  the  final  and  K(r,v)  is  the  initial  electronic  state  potential  energy  as  a  func¬ 
tion  of  the  positions  r-V  =  rj . rv  of  the  .V  nuclei.  Figure  l  illustrates  this  concept  for  a 

diatomic  molecule. 
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We  make  the  approximation  that  the  frequency  of  the  electronic  oscillation  is  so  much 
greater  than  that  of  the  nuclear  motions  that  we  can  evaluate  the  spectral  density  for  clamped 
nuclei,  i.e.  for  r‘v  fixed.  Thus  we  may  write  the  electronic  dipole  moment  at  time  t,  as 

n(t)  -  iifi(ts)  exp[/'n/,(r'v)  t]  (2.10) 

in  which  ji fi  (rv)  gives  the  direction  and  magnitude  of  the  electronic  transition  dipole  moment 
between  states  /  and  /as  a  function  of  the  positions  iN  of  the  nuclei,  i.e.  the  orientation  and 
vibrational  displacement  of  the  molecule,  and  ft  /,  (r'v)  is  the  angular  frequency  of  the  elec¬ 
tronic  transition  dipole  oscillation  at  nuclear  positions  rv  derived  by  Eq.  (2.9)  from  the 
difference  in  potential  energies  of  the  final  and  initial  states  at  positions  tN.  The  only  time  vari¬ 
ation  of  flit)  during  the  evaluation  of  the  spectral  density  is  that  of  exp  (/ft  fit). 

This  can  be  interpreted  as  the  molecular  equivalent  to  the  semiclassical  atomic  "virtual 
harmonic  oscillators”  linked  to  atomic  transitions  which  were  used  by  Landenburg,  Bohr.  Kra¬ 
mers,  Slater,  Bom  and  others  to  describe  the  absorption,  emission,  and  scattering  of  radiation 
by  atoms,  and  which  provided  much  of  the  groundwork  for  the  eventual  development  of  quan¬ 
tum  mechanics.9  The  generalization  to  molecules  is  to  make  the  frequency  of  the  oscillator 
(corresponding  to  the  electronic  transition)  vary  with  the  positions  of  the  nuclei  following  the 
energy  difference  between  the  potential  curves  of  the  final  and  initial  states,  clamping  the  nuclei 
in  position  while  a  spectral  measurement  is  made,  and  then  averaging  over  the  distribution  of 
nuclear  positions  r'v  on  the  initial  potential  surface. 

Thus  evaluating  the  spectral  density10  from  Eqs.  (2.7b)  and  (2.10) 

D[i  -m)  -  12  -M/;(r")l2  5[co  -  n//(r‘v)l  .  (2.1 1) 

in  which  8[w  -  n /j(rv)l  selects  out  the  absorption  frequency  w  equal  to  ft rj( rv). 

Then,  the  cross  section  cr(w)  for  absorption  at  frequency  u  can  be  evaluated  from  Eqs. 
(2.6),  (2.8).  and  (2.11)  as 

<r( at)  —  (  |2  v»/,(rv)l:  8[«  —  ft/,(rv)] )  (2.12) 

in  which  ft ^  is  given  by  Eq.  (2.9),  y  —  (l  —  for  an  equilibrium  system  and  y  —  1  for  a 
nonequilibrium  system  for  which  the  density  of  final  states  pf  is  essentially  zero,  and  ( )  indi¬ 
cates  an  average  over  the  positions  rv  of  the  nuclei  in  the  system.  This  average  may  be 
evaluated  by  molecular  dynamics  for  arbitrary  systems,  and  in  addition  by  Monte  Carlo  or  by 
explicit  integration  over  the  classical  or  quantum  mechanical  distribution  of  nuclear  positions  rv 
for  equilibrium  systems. 

We  illustrate  here  the  calculation  of  electronic  spectra  using  molecular  dynamics  as  we 
wish  to  make  calculations  of  non-equilibrium  transient  spectra.1-2  Thus,  by  computing  the  clas¬ 
sical  trajectories  of  the  nuclei  rv(r)  and  knowing  as  a  function  of  nuclear  positions  the  elec¬ 
tronic  potential  energies  Vf  and  V,  as  well  as  the  electronic  transition  dipole  moment  ii  fl  we 
can  compute  the  electronic  absorption  spectrum  in  the  manner  shown.  We  take  samples  at 
different  times  during  the  trajectories  of  the  nuclei  and  bin  the  instantaneous  power  spectra  into 
appropriate  intervals  in  frequency  to  accumulate  the  averaged  spectrum. 

If  the  system  is  an  isotropic  one,  we  can  average  Eq.  (2.12)  over  all  orientations  to  get4-6 
<r(w)  -  (  |M/l(rv)|2  8h>  -  ft„( rv)J \  .  (2.13) 

C.  Thermal  harmonic  quantum  correction 

A  simple  quantum  correction  can  be  applied  to  the  classical  computation  of  equilibrium 
spectra  to  bring  them  into  closer  agreement  with  quantum  reality.  If  we  make  the  approxima¬ 
tion  that  the  nuclear  motion  of  the  system  can  be  treated  as  harmonic  oscillation,  we  find  that 
the  classical  result  can  be  quantum  corrected  by  simply  scaling  the  temperature  at  which  the 
classical  calculations  are  performed. 
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The  exponential  factor  of  the  diagonal  element  of  the  density  matrix  for  a  harmonic  oscil¬ 
lator  at  temperature  Tis11 


p(x,x:T)  ce  exp 

_  7}  to 


-mm 


x2  tanh  / 


(2.14) 


in  which  /  =  ‘‘  "*■  ,  kB  is  Boltzmann's  constant,  x  is  the  position,  m  is  the  mass,  <u  is  the 

2kg  I 

angular  frequency  of  nuclear  vibration,  H  is  Planck's  constant  divided  by  2ir,  and  tanh  is  the 
hyperbolic  tangent. 

If  we  take  the  correspondence  principle  classical  limit  as  ff—O  of  Eq.  (2.14),  we  find 


!im  exp 

Jl— o 


— -p-x2  tanh  / 


—  exp 


-mwV 


2kBT 


-exp 


FOc) 


k„T 


in  which,  for  a  harmonic  oscillator, 

K(x)  -  ma.V/2. 


(2.15) 


(2.16) 


This  correspondence  principle  result  agrees  with  classical  statistical  mechanics.  Now  we 
compare  the  quantum  and  classical  expressions,  Eq.  (2.14)  and  Eq.  (2.15),  to  calculate  the 
quantum  correction.  Let  T  be  the  physical  temperature  at  which  we  wish  to  compute  the  spec¬ 
tra  and  let  T  be  the  scaled  temperature  for  a  classical  calculation,  designed  to  replicate  the 
quantum  spectra  at  T. 


7*  _  1 
2  ka  tanh  / 


(2.17) 

(2.18) 


Dividing  Eq.  (2.18)  by  T gives  the  quantum  correction  factor  C. 

c  m  r.  _  fo*  i _ _  / 

**  T  2 kBT  tanh  /  tanh  / 


(2.19) 


This  temperature  quantum  correction  is  plotted  in  Fig.  2. 

For  example,  the  scaling  factor  used  to  compute  the  equilibrium  electronic  absorption 
spectrum  for  gas  phase  I2  at  298  K  is  C  —  1.03.  This  means  that  if  we  perform  the  classical 
calculation  at  T  —  307  K  then  we  will  get  the  correct  quantum  spectral  band  contour  in  the 
approximation  of  harmonic  motion  as  if  we  had  done  the  quantum  calculation  at  298  K. 


D.  Symmetry 

For  a  nonequilibrium  system,  the  isotropic  average  of  Eq.  (2.13)  will  not  in  general  be 
correct.  Therefore  we  derive  the  proper  averaging  for  the  important  case  of  pumping  and  prob¬ 
ing  with  linearly  polarized  light.  We  determine  the  average  of  D[i  ’ft]  over  all  orientations  for 
the  anisotropic  and  nonequilibrium  case  in  which  we  prepare  the  sample  by  pumping  it  with 
linearly  polarized  light,  and  then  measure  its  spectrum  by  probing  with  light  polarized  linearly 
in  a  different  direction.  For  a  diatomic  molecule  (symmetric  top)  the  transition  dipole  ft}, 
must  either  lie  along  the  axis  of  the  molecule  or  in  a  plane  at  right  angles  to  the  molecular  axis. 
The  absorption  probability  for  pump  light  is  dependent  on  the  angle  between  the  electric  vector 
i*  of  the  pump  light  and  the  transition  dipole  .ector  ft.,  and  is  greatest  when  (tfl  and  «'  are 
aligned.12  and  thio  the  pump  light  prepares  an  anisotropic  sample  of  excited  molecules,  which 
can  then  be  probed. 

The  angular  distribution  of  excited  molecules  is12 
1(9)  -  -L  [i+0P,(j'  .  *)J 


(2.20) 


f/tanh  f 


O  po  oj  cn 


Figure  2.  Graph  of  harmonic  quantum  correction  factor  T!T  “  //tanh  /.  in  which  T  is  the 
true  temperature  and  T  is  the  quantum  corrected  temperature  at  which  the  classical  calculation 
is  carried  out  and  /  3  fcti2ka  T 


-5- 


where  the  second  order  Legendre  polynomial  is 

Pi(x )  -  3^2~-  (2.21) 

and  9  is  the  angle  between  the  electric  unit  vector  «'  of  the  pump  light  and  the  unit  vector  h 
along  the  axis  of  I2.  The  (3  term  is  referred  to  as  the  asymmetry  parameter.  When  the  transi¬ 
tion  dipole  vector  it/i  points  along  the  molecular  axis,  /3  -  2,  a  cosine-squared  distribution. 
When  /if i  is  perpendicular  to  the  molecular  axis,  -  —1,  a  sine-squared  distribution. 

We  would  like  to  compute  (  Dli'/i] }  which  determines  the  cross  section  for  absorption 
by  Eqs.  (2.6)  and  (2.8)  in  which  ( )  indicates  an  average  over  all  orientations  for  the  distribu¬ 
tion  of  the  molecules  which  have  been  excited  by  the  pump  light.  We  can  write  this,  using  the 
notation  of  Eq.  (2.7),  as 


(  0(i •*»!)*(  Dlfi^iij] )  (2.22a) 

ij 

-I«,  (2.22b) 

<J 

”  X  A,j  (2.22c) 

ij 

in  which  the  subscripts  /  and  prefer  to  the  Cartesian  components  of  i  and  /*(/  '•  A,,  is  the 

average  of  over  the  orientational  distribution  of  molecules  excited  by  . .  np  light. 

We  want 

Ay  -  (  )  (2.23a) 

-JV*  1(9)  (  ),  (2.23b) 

-  ‘  A)I '  Db,#j\  \  (2.23c) 


in  which  h  is  the  unit  vector  along  the  axis  of  I2  at  the  time  of  excitation.  1(9)  is  the  probabil¬ 
ity  distribution  for  excitation  and  ..  D[n,,jij]  )  is  the  average  conditional  on  the  axis  of  the 
molecules  being  fixed  in  a  certain  direction  h  when  the  pump  light  is  absorbed. 

Using  rotational  symmetry  (see  Appendix  A) 

,  -  Ai0)  6,  +  A'»  (2.24) 

so  that  for  molecules  with  intemuclear  axes  aligned  along  the  :  axis 

\  -  A,m  +  A,2>  (2.25) 

DU,#,  1  f  -  DU,#, I  }  -  A,0)  -  y  A,2>  (2.26) 


and 


y  L  1 )  -  A 

*  i  - 


(0) 


Combining  Eq.  (2.24)  with  the  following  (see  Appendix  B) 

f  d1  n  .  .  1  . 

J  4t  “  3  5,1 

f  4“  ",  n,  hk  h,  -  - jj  [5,;5W  +  5*8,,  +■  M,*)  . 


Art  •”  ">  3  " 

1  d1  h  . 

4-rr  15 

Eq.  < 2.23c)  can  be  evaluated  as 


(2.27) 

(2.28) 
(2.29) 


A.,  -  .4'°'  +0  A' 


10 


(2.30) 
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Finally,  we  write  our  average  from  Eq.  (2.22c)  as 


(  DU  -fi  1  }  -  Z  e,  ij 


i<0>  s  -l.  a  3  j  1  *  >  S-V 


6,,+PA' 


10 


-  Ai0>  +  0  A 


m  3(«*  •  «)2  -  1 
10 


(2.31a) 

(2.31b) 


-A{0)  +  &■  A(:>  •«) 


(2.31c) 


in  which  i'  is  the  electric  unit  vector  of  the  pump  light  and  «  is  the  electric  unit  vector  of  the 
probe  light. 

Therefore,  to  compute  the  absorption  spectrum  for  probe  light  polarized  along  I  by  a  sys¬ 
tem  prepared  by  pump  light  polarized  along  5' , 

<r(w)  -  lAm  +  4  A,2>  P2(i'  •  €l|  (2.32) 

Tien  5 

in  which 

A‘0)  -  j  (  |/t/,(rv)l2  SU  -  n/,(rv)l )  (2.33) 

A (2)  -  (  [a  •#»/,(rv)lJ  8[«  -  n/((rv)J },  -  ^(0'  (2.34) 

where  /3  is  the  asymmetry  parameter  defined  above  for  the  pump  light,  ft/,(ry)  is  the  transition 
dipole  vector  at  nuclear  positions  r  v  for  the  probe  light,  and  {  [«  /,  (rv)l2  .  is  the  average  of 

the  h  component  squared  of  the  probe  transition  moment  vector  at  nuclear  positions  rv  for 
molecules  with  symmetry  axis  aligned  along  h  when  the  pump  light  is  absorbed. 


III.  RELATION  TO  OTHER  TECHNIQUES 

The  almost  classical  technique  we  illustrate  here  is  related  to  other  sermclassical  tech¬ 
niques  for  computing  electronic  absorption  spectra.7-3  This  approach  is  related  to  a  short  time 
limit  of  the  wave  packet  time  dependent  method.3- 13-16  Thus,  it  is  appropriate  for  liquid  solu¬ 
tions  with  rapid  randomization  (which  can  also  be  treated  by  wavepacket  techniques17  ).  for 
excitation  to  quickly  dissociating  states,  and  for  the  band  contours  of  systems  with  vibrational 
and  rotational  structure,  the  details  of  which  would  be  given  by  longer  time  wave  packet  evolu¬ 
tion. 

In  this  section  we  show  the  equivalence  of  the  binning  method  used  for  calculations  in 
this  paper  and  the  semiclassical  methods  of  Lax,18  and  Tellinghuisen  and  Moeller19  and  Lee.30 
Lax's  discussion18  of  the  relation  of  "exact",  and  classical  Franck-Condon  spectra  applies  equally 
to  the  comparison  of  "exact"  wavepacket  and  molecular  dynamic  binning  spectra,  and  his  semi¬ 
classical  Franck-Condon  spectra  are  approximated  by  our  harmonic  quantum  correction.  The 
Landau  -  Zener  -  Stuckelberg  surface  crossing  probabiiny.  in  a  form  similar  to  the  Tully  -  Pres¬ 
ton  surface  hopping  method,21  is  shown  to  give  the  same  results  as  the  other  methods  when 
the  lower  surface  is  raised  by  Tko.  giving  rise  to  intersection  is)  of  the  ground  and  excited  state 
surfaces,  which  can  be  treated  as  surface  crossings  with  a  strength  given  by  the  usual  dipole 
coupling.  The  binning  method  given  in  this  paper  is  considerably  simpler,  both  conceptually  and 
computationally,  than  the  other  methods  mentioned  above,  so  it  is  important  to  establish  the 
equivalence.  In  addition,  the  extension  of  the  molecular  dynamic  binning  method  to  nonequili- 
bnum  systems1-2  is  particularly  clear  and  simple. 


A.  The  molecular  dynamic  binning  method 

Consider  Fig.  3  which  shows,  for  a  two  dimensional  case,  a  contour  plot  of 
Vf(ts)  —  V,  (rv),  together  with  a  particular  trajectory  on  the  V,  potential  surface.  This  trajec¬ 
tory  is  plotted  as  dots  spaced  at  equal  times,  which  are  the  sample  times  used  to  inquire  as  to 
which  bin  the  trajectory  is  in.  The  space  between  a  pair  of  contour  lines  corresponds  to  a  fre¬ 
quency  range 

a»*  —  ^  <u  ^  <uk  +  ,  (3.1) 

where  otk  —  k  8<o,  i.e.  the  ktH  bin.  (The  contour  intervals  are  then  8E  —  Ti  6o>).  In  Fig.  4  we 
magnify  a  portion  of  Fig.  3  and  show  a  few  more  details.  The  amount  of  time  spent  between 
the  two  contour  tines  shown  is  proportional  to  the  contribution  this  particular  trajectory  seg¬ 
ment  will  make  to  the  corresponding  spectral  frequency  bin.  The  time  r  spent  between  con¬ 
tours  is  the  distance  D  traversed  divided  by  the  velocity,  or 


D 

(p/m) 


md 
p  sind 


in  which  p  is  the  magnitude  of  the  momentum,  pL  —  p  sind  is  the  magnitude  of  the  momentum 
perpendicular  to  the  contour  lines,  and  m  the  mass.  The  perpendicular  separation  d  between 
the  contours  is 


where 


UFJ  -  |V[  F/(rv)  -  K(rv)]| 


evaluated  in  the  vicinity  of  the  crossing  region.  is  the  force  change  between  the  final  and 
initial  potential  surfaces  perpendicular  to  the  contour  lines  of  V,  -  V,  at  the  intersection 
V,  +■  -  Vf.  Taking  the  rate  of  crossing  per  unit  time  from  the  initial  to  the  final  surface 

to  be 

*  -  .  (3.5) 

7r  5  to 

in  which  W  is  a  coupling  parameter,  the  total  probability  P  of  crossing  from  the  initial  to  the 
final  surface  while  on  the  trajectory  in  the  region  of  the  bin  is 


P-kr - 


2itW2  1  I  mf\  8(o 


8<u  Pi  UF. 


2k  W2 

n 


Pi  UFJ 


(3.6b) 


Nothing  essential  depends  on  our  assumption  of  two  degrees  of  freedom,  and  Eq.  (3.6)  is  gen¬ 
eral. 


B.  Equivalence  with  Landau  •  Zener  -  Stuckelberg,  Tully  •  Preston  (LZSTP)  method 

Eq.  (3.6)  is  just  the  Landau  -  Zener  -  Stuckelberg  rate  for  radiationless  transtion  surface 
crossing21- 22  jn  the  case  0f  diabatic  surfaces  which  intersect  with  small  coupling  W.  (The  adia¬ 
batic  crossing  probability  would  be  1  -  P.)  The  multidimensional  picture  of  a  hopping  trajectory 
is  due  to  Tully  and  Preston.21  The  calculation  of  the  rate  for  crossing,  as  opposed  to  following 
the  trajectones  after  they  have  crossed  (which  we  do  not  need  to  do),  has  been  recently  carried 
out  in  the  context  of  radiationless  transitions  (cu  -  0)  by  Heller  and  Brown.22  (who  did  not  use 
the  simplified  binning  method)  and  the  results  are  excellent. 

For  coupling  IF caused  by  dipole  radiation,  we  have,  for  unit  field  strength. 


Figure  3.  A  two  dimensional  contour  plot  of  V,  (r  v)  —  V,  (rv)  showing  a  particular  trajectory 
on  the  V  potential  surface  plotted  as  solid  circles  spaced  at  equal  times.  The  mset  box  (dashed 
lines)  is  magnified  in  Fig.  4 


I 


Figure  4.  A  magnified  view  of  the  inset  box  of  Fig.  3.  0  is  the  distance  traversed  between  adja^ 
cent  contour  lines  and  d  is  the  perpendicular  distance  between  the  contour  lines. 
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H/2(rv)  -  i  |M/l(rv)|J ,  (3.7) 

ft 

in  which  n  is  the  index  of  refraction. 

C.  Equivalence  with  "classical"  method 

It  is  possible  to  arrive  in  several  different  ways  at  a  semiclassical  expression  for  the 
absorption  cross  section  which  involves  only  classical  (possibly  quantum  corrected)  phase  space 
densities.  Lax18  was  apparently  the  first  to  derive  it.  Lee20  calls  it  the  "semiclassical"  cross  sec¬ 
tion,  and  Tellinghuisen  and  Moeller19  call  it,  as  we  do  here,  the  classical  absorption  spectrum. 
Indeed,  we  arrived  at  their  formula  in  Section  II  in  yet  another  way,  i.e.  the  classical  oscillating 
dipole  analogy  of  the  old  quantum  theory. 

To  show  the  equivalence  of  Eq.  (2.13)  and  Eq.  (25)  of  Tellinghuisen  and  Moeller19  to  the 
binning  method  and  LZSTP.  we  need  to  relate  the  individual  trajectory  picture  used  in  Sects.  Ill 
A  and  B  to  the  coordinate  space  integral  expression,  Eq.  (2.13). 

The  concept  of  the  binning  method  is  easiest  to  work  with.  We  imagine  pointwise  sam¬ 
pling  the  appropriate  phase  space  density  p(rv,pv),  which  may  be  equilibrium  or  nonequili¬ 
brium.  We  merely  check  to  see  what  bins  the  phase  space  points  are  in,  and  assign  them 
appropriately,  thus  building  a  histogram  of  the  absorption  profile.  There  is  no  need  to  run  the 
trajectories,  other  than  as  a  means  to  compute  p( rv,pv),  since  in  the  short  time  required  to 
cross  a  bin  the  overall  phase  space  density  is  not  changing  much  (or  for  a  stationary  p(rv.p  v) 
not  changing  at  all).  Also,  momentum  is  irrelevant  as  to  which  frequency  bin  the  phase  space 
point  belongs.  Thus,  the  total  rate  of  crossing  induced  by  the  electromagnetic  field  is  the  bin 
count  over  all  of  phase  space,  or 

k  -  f  dr*  dp  v  p(r  v,p  v)  lK2(rv)  Alw  -  n„]  (3.8a) 

37f  ooi  J 

-  /  dxs  P( rv)  lM/,(rv)lz  A(cu  -  fl/()  .  (3  8b) 

where  P( rv)  is  the  probability  as  a  function  of  position,  the  factor  of  1/3  is  introduced  for 
spherically  averaged  molecular  orientations,  and 

j  1  ~'/i8w  <  X  < 

A(T)  -  j  q  otherwise. 

Now  take  the  limit  as  5<u— 0 

k  "  f  dt'S  PiT‘W)  ^/'(r'V)|J  8(w  “  n/-}  •  (3  8c) 

To  relate  this  rate  to  a  cross  section,  we  note  that  for  unit  field  strength  the  photon  flux  is 
c/2ir<u.  and  cr  •  flux  —  rate,  thus 

<T<o,)  -  f  dry  P( rv)  lM/,(rv)lJ  5(w  -  Q/()  <3.8d) 

-  4E2L  |M/,(rv)|2  5(a,  -  fl/#(rvH)  •  (3.8e) 

Eq.  (3.8d)  is  the  multidimensional  form  of  the  diatomic  work  of  Tellinghuisen  and  Moeller,19 
<see  also  Lax18  and  Lee20  )  and  Eq.  (3  8e)  is  Eq.  (2.13)  of  this  paper.  As  this  equation  is 
derived  from  the  binning  method,  the  equivalence  is  thus  established  between  LZSTP.  the  bin¬ 
ning  method  and  the  "classical  spectrum". 
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0.  Relation  to  reflection  method 

It  is  natural  to  seek  the  relation  to  the  reflection  approximation,  probably  the  most  used 
technique  for  estimation  of  absorption  spectra  in  diatomic  molecules.  Heller14  and  Lee.  Brown 
and  Heller23  have  given  the  correct  polyatomic  generalization  of  the  reflection  method.  It  is 
related  but  not  equivalent  to  the  methods  of  Sects,  ill  A  •  C.  For  low  lying  vibrational  levels  of 
V,  and  steep  upper  surfaces  the  reflection  method  should  be  superior  to  classical  absorption 
techniques. 

As  Tellinghuisen  and  Moeller19  have  noted,  the  reflection  approximation  has  its  draw¬ 
backs,  especially  for  moderate  and  high  temperatures,  since  each  quantum  wavefunction  with 
significant  Boltzmann  populations  must  be  known  and  used.  The  higher  vibrational  states  give 
rise  to  inaccurate  results  (for  a  discussion  of  this  see  Heller14  ).  and  the  method  is  therefore 
suspect  at  high  temperatures,  where  it  is  also  most  difficult  to  apply.  Here,  the  classical  spec¬ 
trum  methods  are  distinctly  advantageous. 

The  reflection  approximation,  in  the  present  notation,  is 

<r( «)  -  ^£2.  I /  dts  |'Mr")|2  lM//(rv)l2«(w  -  Wf(t»)  -  £„))  .  (3.9) 

The  main  distinction  between  Eq.  (3.9)  and  Eq.  (3.8)  is  the  replacement  of  K,(rv)  by  the 
energy  of  the  /*  vibronic  state  on  the  initial  surface  4  and  the  enumeration  of  each  initial  wave 
function  yY,j(is). 

E.  Wavepacket  methods 

Wavepacket  dynamics  can  also  be  used  to  compute  electronic  spectra3- 17  and  provides  a 
route  to  more  exact  solutions.  In  a  semidassical  version  the  method  is  efficient  but  not  as  fast 
as  the  more  classical  methods  (especially  the  binning  method)  discussed  here.  However,  use  of 
the  coherent  wavepacket  approach  permits  accurate  rendition  of  the  spectrum,  including 
interference  oscillations  and  spectral  structure  not  available  in  any  of  the  methods  mentioned 
above.  There  are  no  restrictions  to  low  or  high  temperatures,  but  for  highly  anharmonic  poten¬ 
tials  approximate  wavepacket  propagation  methods3  can  have  significant  errors  for  higher  reso¬ 
lution  spectra  and  more  exact  propagation  may  be  required. 13 

IV.  INPUTS  TO  Ij  CALCULATION 

In  this  section,  we  test  the  binning  molecular  dynamics  method,  computing  the  band  con¬ 
tour  of  the  visible  absorption  spectrum  of  room  temperature  gas  phase  Ij  molecules,  and  com¬ 
pare  the  results  to  experimental  measurements. 

A.  Potential  Energy 

Figure  5  shows  the  potential  energy  curves  for  the  iodine  ground  X  0*('X)  state  and  the 
excited  A  ltt(Jn).  BO^On)  and  B"lM('n)  states.  The  ground  X  state  is  taken  from  a  semi- 
empirical  potential  due  to  Matzen.  Calder  and  Hoffman24  and  from  the  experimental  RKR 
values  of  Coxon.25  The  A  state  and  the  B"  repulsive  state  are  from  Tellinghuisen,26-27  and  the 
B  state  is  from  Barrow  and  Yee28  and  Luc.29 

For  the  ground  state  for  imemuclear  distance  R  <  0.323  nm  and  for  the  entire  B"  state 
the  potential  energies  are  easily  calculated  from  potential  functions  V(R)  cited  above.  The 
ground  state  for  R  >  0.323  nm,  the  A  state  and  the  B  state  potential  energy  values  are  gen¬ 
erated  by  the  Aitken-Neviile  technique  30  of  iterated  linear  interpolation  of  the  RKR  turning 

points,  a  set  of  n  points  V<R,),  i  -  0 . n  where  V(R,)  is  the  potential  energy  at  intemudear 

distance  R,.  Interpolation  methods  use  function  points  to  determine  a  unique  polynomial 
P„  ( R )  which  satisfies  the  constraints 

P„(R,)  -  V(R,),  (-  0 . n.  (4,1) 

The  Aitken-Neviile  method  uses  repeated  linear  interpolation  to  compute  the  value  of  Pn(X) 


INTERNUCLEAR  DISTANCE  (nm) 

Figure  5.  The  potential  energy  curves  for  the  ground  X  0p'£'.  A  IJ-’FP-  B  0~<JrP-  jnd 
states  used  in  the  calculation  of  the  electronic  absorption  spectrum  of  iodine.  The 
dots  show  the  RKR  turning  points  for  every  fifth  vibrational  level. 
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for  X in  (R0,R„)  subject  to  the  constraints  that  P„(R,)  -  V(R,).  Here  AMs  an  intemudear  dis¬ 
tance  for  which  a  value  of  the  potential  energy  V(X)  is  required.  Ten  function  points 

K(/?,),  /  —  m-9 . m,  are  used  in  the  interpolation  to  compute  P„(X)  =  V(X).  These  points 

are  a  subset  of  the  complete  set  of  RKR  points.  The  range  (/?„,_$, Rm)  is  chosen  such  that  X 
lies  as  close  as  possible  to  the  median  value.  In  the  case  where  AMs  close  to  R0  or  R„,  where  n 
is  the  total  number  of  RKR  points,  there  are  still  ten  points  used  in  the  interpolation. 

The  upper  portions  of  the  inner  repulsive  branch  of  the  A  state  and  the  B  state  are 
approximated  by  potentials  of  the  form  V(R)  -  A+C/R”.  In  the  case  of  the  B  state  repulsive 
branch  C  -  8.2805  x  10-4,  n  «■  12  and  A  —  164.60.  The  smooth  curve  is  attached  at  v’  -  34. 
For  the  A  state  repulsive  branch  C  -  3.538  x  10~5,  n  —  11,  A  —  104.87  and  the  curve  is 
attached  at  v'  -  35.  Here  the  units  of  A  and  C  are  such  that  V(R)  is  in  kJ  mole-1  when  R  is  in 
nanometers. 


B.  Transition  dipole 

The  transitions  strengths  Im/jI2  for  the  A— X  transition  and  the  B"—X  transition  are 
taken  to  be  constant  as  a  function  of  intemudear  distance,  with  values27  lp/,1*  x  “  0.453 
x  lO-^CW  (0.0407  D2)  and  l*/,l$._x  -  1.5135  x  10'60  C2m2  (0.136  D2).  in  which  D  is 
Debyes  (eA).  For  the  B— X  transition  the  electronic  transition  moment  versus  R  cen¬ 

troid  is  derived  over  the  range  R  “  0.25  to  R  —  0.29  nm  from  the  transition  dipole  data.27  The 
data  points  for  Im/,1  from  Tellinghuisen  27  are  used  to  construct  two  linear  segments  for  the 
intervals  R  *■  0.25  to  0.2685  nm  and  R  *  0.2685  to  0.29  nm.  The  resulting  Ijt/,  I  versus  inter- 
nudear  distance  functions  for  the  three  transitions  are  shown  in  Fig.  6. 


V.  COMPUTATIONAL  METHOD 

To  compute  the  classical  trajectories31  rv(r)  -  rj(r) . r.v(r)  we  first  choose  an  initial  set 

of  coordinates  and  momenta  for  the  N  atoms  in  our  system  and  then  describe  the  atomic 
motion  by  integrating  Newton’s  second  Jaw 


av 

dr, 


-  F, 


m, 


fh. 

dt 2 


i  “  1  , 


(3.2) 


in  which  V  -  V(rv)  is  the  potential  energy  of  the  atoms  at  positions  r,,  .  .  .  ,r  v,  F,  «■  F,(rv)  is 
the  force  on  the  ith  atom,  and  m,  is  the  mass  of  the  /  th  atom.  The  time  step  size  is  1.0  femp- 
tosecond.  A  modified  Verlet  integration  algorithm  is  used32'34  along  with  minimum  image 
periodic  truncated  octahedral  boundary  conditions35  to  reduce  edge  effects. 

Before  we  actually  make  spectral  calculations,  we  equilibrate  the  system  at  the  desired 
temperature  7*by  integrating  forward  in  time,  stopping  at  intervals  to  choose  a  new  set  of  velo- 
dties  selected  at  random  from  a  Maxwell-Boitzmann  distribution  consistent  with  the  tempera¬ 
ture  T.  Then,  to  sample  more  of  configuration  space,  we  again  pick  new  sets  of  velocities 
between  sets  of  trajectories  used  for  the  spectral  calculations,  and  average  the  resulting  spectra. 

The  classical  trajectories  for  a  system  of  100  non-interacting  I2  molecules  were  computed 
for  a  time  period  of  6  nanoseconds. 


VI.  Ij  EQUILIBRIUM  GAS  PHASE  SPECTRUM 

Figure  7  shows  the  equilibrium  gas  phase  spectrum  of  I2  as  calculated  by  molecular 
dynamics  and  the  gas  phase  experimental  data  from  Tellinghuisen27  superimposed  as  solid  cir¬ 
cles.  The  solid  line  is  the  quantum  corrected  gas  phase  spectrum  while  the  dashed  line  is  the 
uncorrected  gas  phase  spectrum.  There  is  very  good  agreement  between  the  theoretical  and 
experimental  electronic  absorption  band  contours  of  I;  and  the  quantum  correction  is  seen  to  be 
small  for  this  case  of  low  ground  state  vibrational  frequency. 
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figure  6.  Transition  dipole  moment  magnitudes  >.,(/?>!  as  functions  of  1:  intemuclear  dis¬ 
tance  R.  For  the  A— X  and  the  B"— X  transitions  ihe  function  is  taken  to  be  a  constant  while 
for  the  B— X  transition  it  is  shown  as  a  function  of  R  centroid.  "Tie  solid  circles  are  the  data 
points  from  Teilinghuisen. 
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VII.  DISCUSSION  AND  CONCLUSION 

We  have  shown  that  the  band  contour  of  the  electronic  absorption  spectrum  of  a  gas 
phase  equilibrium  system  can  be  computed  by  a  nearly  classical  mechanical  approach  and  have 
used  I;  as  an  illustration.  Utilizing  potential  energy  and  transition  moment  functions  from  the 
literature,  we  have  been  able  to  calculate  the  electronic  absorption  spectral  band  contour  to  a 
remarkable  accuracy.  Furthermore  we  have  shown  that  a  quantum  correction  can  be  applied  to 
the  classically  computed  spectrum,  but  that  for  the  massive,  weakly  bound  I*  molecule  this 
correction  is  quite  small.  Since  classical  molecular  dynamics  are  the  basis  of  these  calculations 
one  can  have  confidence,  given  the  agreement  of  computed  and  measured  spectral  contours, 
that  the  classical  nuclear  motions  as  simulated  are  a  reasonable  classical  representation  of  the 
actual  nuclear  time  evolution  of  the  real  system. 

Another  very  different  way  to  compute  such  spectra  ab  initio  from  molecular  dynamics 
would  be  to  use  a  quantum  force  classical  trajectory  (QFCT)  method36  in  which  we  follow  a 
classical  trajectory,  computing  at  each  time  step  the  electronic  transition  dipole  vector  p /,  (r  v) 
and  the  potential  energies  V/(ts)  and  K,  (rv)  of  the  final  and  initial  states  in  order  to  calculate 
the  transition  dipole  frequency  ft /,  (rv),  and  to  bin  the  results  to  accumulate  the  spectrum. 

Calculations  of  electronic  band  contours  with  the  method  presented  here  are  particularly 
appropriate  for  the  liquid  state  in  which  the  rotational  and  much  of  the  vibrational  structure  is 
washed  out  by  dephasing  due  to  averaging  over  solvent  molecule  initial  conditions  and  by 
different  interactions  with  solvent  molecules  in  the  ground  and  excited  states.17  We  have  illus¬ 
trated  liquid  state  calculations  elsewhere1-2  using  the  technique  presented  in  this  paper.  Due  to 
the  simplicity  and  close  connection  to  classical  trajectories,  we  can  also  easily  carry  out  non- 
equilibrium  time  dependent  calculations  as  we  have  already  shown1-2  to  compute  the  transient 
absorption  during  chemical  reaction?  in  solution. 
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APPENDIX  A:  PROOF  OF  EQ.  (2.24) 

We  wish  to  explore  the  consequences  of  rotational  symmetry  for  the  structure  of  the 
quantity 

A,,(h)  —  I  n  ■  (Al) 

For  any  (proper  or  improper)  rotation  R  -  [/?,,].  rotational  invariance  implies  that 

A„(R  *)  -  I  Mm  Akl(H).  •  <A2) 

w 

as  can  be  seen  by  recognizing  that  for  any  pair  of  vectors  a  and  b,  Y.A„lh)a,b,  is  invariant 
under  rigid  rotation  of  the  triple  a.  b,  h. 

First  consider  the  case  h  -  2  and  R  a  rotation  by  <b  about  the  :  axis.  Differentiating  Eq 
(A2)  with  respect  to  <b  and  then  setting  <b  to  zero  yields 


^12(J)+4ji<2)  A  22<2)  — ,42jf2) 

—  ,-f  1 2  ( i  ) — .4’1 12)  —^13(2) 


A  32(2)  ~A]  [  ( 2 )  0 


-  0. 


(A3) 


from  which  we  conclude  that 
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Am-  KAi2)  -A'"  0 

A  (2)-  A{"  A(0)-'/2Ai2)  0  (A4) 

0  0  A,0)  +  A,2> 

for  appropriate  constants  A(‘\  i  -  1,2,3. 

Next,  invoke  Eq.  (A2)  for  k  -  2  and  R  the  improper  rotation  (reflection  in  the  xz  plane). 
We  are  not  assuming  that  the  excitable  molecule  has  an  inversion  center,  any  more  than  we 
assumed  the  excitable  molecule  to  be  spherically  symmetric  in  the  preceding  paragraph. 

1  0  0 

R  -  0-1  0  .  (A5) 

0  0  1 

Using  Eqs.  (A5)  and  (A3)  in  (A2)  yields  An)  -  —AiU,  so  A11’  -  0. 

We  can  write  Eq.  (A3),  with  A(1)  •  0,  in  the  form 

A,j(2)  -  A,0%  +  Am—  !^±L.  (A6) 

Now  choose  any  rotation  R  such  that  R  2  •  h,  i.e.,  R,j  -  hi%  i  -  1,2.3.  Using  Eqs.  (A2)  and 
(A6),  we  have 

A,j(i I)  -  A0(R  2) 

-  I R*Rj,Au(i) 

u 

-  +  a™ 

Id  * 

-  A'%  +  A(2>  rnrnC*'L ,  (A7) 

where,  in  the  last  step  we  used  the  orthogonality  of  R 

I*, **,/«« -I*, */?,*-«*  (A8) 

and  R  2  -  k  ( so  that  -  »,•).  With  the  definition  Eq.  (Al),  Eq.  (A7)  is  identi¬ 

cal  to  Eq.  (2.24) ,  which  we  set  out  to  prove. 

APPENDIX  B:  PROOF  OF  EQS.  (2.28)  AND  (2.29) 

Define  the  linear  map 

\f:  £3®£3  —  Es.  (Bl) 

where  E"  is  the  real  linear  space  of  ^-tuples,  by  the  matrix 


.V/  is  the  matrix  that  couples  two  p-states  to  a  Estate.  From  (B2)  it  is  easy  to  read  off  that 
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MMr  —  1,  (the  identity  matrix  on  f5)  and  that  P  -  Mr M  is  the  projection  (on  d-states) 
whose  matrix  is 

P.JM  ”  to  —  yS,y5t/).  (B3) 

The  point  of  introducing  M is  that,  for  any  rotation  R  there  is  a  5  x  5  matrix  H/t)  such  that 
M  R®R  -  T{R)  M,  (B4) 

and  T(A)  is  a  unitary  irreducible  representation  of  the  rotation  group. 

1.  The  existence  of  T(/?)  and  its  being  an  irreducible  representation  follow  from  the  obser¬ 
vation  that  the  polynomials 

(BS) 

•j 

are  Estate  wave-functions. 

2.  The  unitarity  of  T(R)  follows  from 

T(R)T(R)r  -  M  R®R  Mt  M  /T1®/?”1  Mt 


«  M  /?«/?  P  R-'QR-'  Mt  (B6) 

-  M  P  Mr  -  MMt  MMt  -  1, 
thanks  to  the  consequence 

R9R  P  /r'®/rl  -  P  (B7) 

of  the  orthogonality  of  R. 

We  are  now  prepared  to  evaluate  rotational  averages  (2.28)  and  (2.29).  Using  (B3). 

MTM  A®A  -  P  A®a  -  A®a  -  -~8.  (B8) 

so,  with  the  definition 

v(A)  -  M  *®#,  (B9) 

and  the  observation,  which  follows  from  (B4)  and  (B9). 

v(r  a )  -  m  r®r  a®a  -  n* )  .v/  a®a  -  n/t>  v(a>.  (bio) 

we  have  a  decomposition  of  a®  a 

a®a  -  y«  +  mt  v(a>  (bid 

into  components  transforming  by  the  one  and  five  dimensional  irreducible  representations  of 


the  rotation  group.  The  strategy  for  evaluating  averages  /(A)'  over  the  unit  sphere  is  to 
replace  them  with  averages  over  the  rotation  group  f(R  z),  and  use  the  orthogonality 
theorem  of  group  representation  theory  T7  In  this  way,  we  find 

v<a»-  n/?) ;  v(f)-o,  <Bi2) 

which,  with  (Bll),  proves  (2.28),  and 

v„(k)va(k))-Z  raJR)raJR)  vji)  vjz) 

*** 

<  B 1 3 ) 

m 1  —  j5 


which,  with  <B11).  yields 
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(My***/)  -  l  8 „hu  +  'Z.Ma.IJMB.ki  (vjh)va(h)) 

y  at 

+  (Bl4) 
a 

”  ^  8/;5*/  +  P,j.kl- 

When  (B14)  is  combined  with  (B3),  we  obtain  (2.29). 
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